// 1D diffusion equation
#ifndef DIFFUSION_EQU_1D_GAS_PHASE_H
#define DIFFUSION_EQU_1D_GAS_PHASE_H

#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;
#include "../../../ThirdPartyCode/Numerical Recipes in CPP 2nd Ed/0_Matrix_and_Vector/nrtypes.h"

#define LAYER_PER_SEC_TO_MOLE_PER_SEC (1/5.7423e4)
#define MOLE_PER_SEC_TO_LAYER_PER_SEC (5.7423e4)

class diff_equ_1D_gas_phase_rec
{
	// Model settings and initial values
	double	Z_MAX;
	double	P;
	double  K;// The reaction model
	Vec_DP	*vecT;
	double	X_ini;		// Assume constant initial concentration
	double	Xin_ini;	// Initial values of inlet concentration

	// Simulation settings
	DP		dz;
	double	t_begin;	
	double	t_end;
	double	dt_rec;
	DP		dt;
	char*	OutputFileName;
	ofstream fout;
	ofstream fout_detail;
	char	FileName2[256];
	int		Detail;
	// I/O
	double	x_BC1;

	// Internal variables
	int		znode;
	Vec_DP	*zARR;
	int		t_node;
	int		t_idx;
	double	*t_arr;
	int		Repeated;
	
	Vec_DP	*x;		// Molar concentration
	double	*Xs_sum;
	double	*Xs2_sum;
	double	*DepRate_sum;
	double	*DepRate2_sum;
	
	Mat_DP *A;
	Mat_DP *LU;
	Vec_INT *indx;
	DP		d;
	Vec_DP *b;

	double diffusivity(double T,double P);
public:
	DP		time;
	double	DepRate;
	double	Xs;
	// 2009.08.01 Gangshi: Add iDetail to regulate the output, 0: no detail files; 1: detail files
	diff_equ_1D_gas_phase_rec(double iK,
					 double iZ_MAX,
					 double iP,
					 double iT_reactor,
					 double iT_surface,
					 double iX_BC1,
					 double iX_IC,
					 double it_end = 10,
 					 double idz = 0.001,
					 double idt = 0.01,
					 double	idt_rec = 0.1,
					 char*  iFileName = "diff_equ_1D_gas_phase_rec.dat",
					 int iDetail = 0);
	~diff_equ_1D_gas_phase_rec();
	double update(double Xin);
	int	   notfinish();
	void get();
	void reset();
	double getXs();
	double getXi();
	double getDepRate();
	void setXi(double iXi);
};

void ludcmp(Mat_IO_DP &a, Vec_O_INT &indx, DP &d);
void lubksb(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b);

#endif